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Abstract. We discuss a fast cross-Wigner transform based technique for 
detecting gravitational wave bursts, and estimating the direction of arrival, using a 
network of (three) non co-located interferometric detectors. The performances of 
the detector as a function of signal strength and source location, and the accuracy 
of the direction of arrival estimation are investigated by numerical simulations. 
The robustness of the method against instrumental glitches is illustrated. 



1. Introduction 

The next generation of interferometric detectors, of gravitational waves (hencefortlr 
GW) including AdLIGO [J, AdVirgo and GEO-HF % hopefully to be followed 
soon by LCGT [1] and ACIGA [5], and eventually by ET [6j, is expected to observe 
tens of events per year, opening the way to gravitational wave astronomy [7]. 
Identifying the direction of arrival (henceforth DOA) of the signals, and retrieving 
their shapes, will be a primary task in reconstructing the physics of the sources and 
their environments. 

The possibility of retrieving the DOA from independent estimates of the signal arrival 
time at each detector was first suggested in j8j, and further discussed in Saulson seminal 
book It was shown that three- interferometers are sufficient to retrieve the DOA up 
to a mirror-image ambiguity which can be solved in principle from knowledge of the 
detectors' directional responses. This method, often referred to as triangulation was 
further elaborated by Sylvestre [lU], Cavalier et al. [TT], and Merkovitz et al. [H]. In 
[11] a Gaussian distribution was assumed for the (independent) arrival time estimation 
errors, and a minimization algorithm was accordingly proposed for retrieving the 
DOA, in the maximum likelihood spirit. In [12] it was shown that this method is 
affected by a systematic bias in the estimated DOA, a possible technique for removing 
the bias was discussed, and ampiitude consistency tests for removing the mirror-image 
ambiguity were suggested. Fairhurst developed a similar analysis of the effect of arrival 
time estimation errors on the DOA estimation accuracy, for the special case of chirping 
signals, including waveform and calibration errors |13j . |14j . 

DOA estimation algorithms are already implemented in the coherent LIGO- Virgo 
pipelines for GW burst (henceforth GWB) detection [TS], [T^]. DOA estimation in 
coherent network data analysis, was studied first by Krolak and Jaranowski [T7], and 
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then by Pai et al. [18j , as part of the waveform parameter estimation problem, with 
specific reference to chirping waveforms from coalescing binaries, in a Gaussian noise 
background. The conceptual foundations of coherent data analysis for unmodeled 
waveforms were laid out by Flanagan and Hughes |19) . and further developed by 
Klimenko et al. [10]- EB- Giirsel and Tinto [32] first suggested the possibility of 
retrieving the DOA for un modeled signals using null-streams. This concept was 
analyzed in depth by Schutz and Wen [23], and further exploited by Chatterji et 
al. [23]. A Fisher-matrix based analysis of arrival time estimation error in coherent 
network detection of modeled as well as unmodeled signals was made in Wen et al. 

m- 

In this paper we capitalize on the time-shift and localization properties of the cross- 
Wigner-Ville (henceforth XWV) transform to introduce a new and conceptually simple 
GWB detection and DOA reconstruction algorithm, using a network of non co-located 
interferometric detectors. 

The Wigner-Ville transform is a well known powerful tool for the analysis of non- 
stationary signals [SS], whose potential in GW data analysis, has been highlighted 
by several Authors, under different perspectives [22]- [29]. Here we suggest its 
possible use as an effective tool for detecting GWBs, and estimating their DOA, 
which offers nice features in terms of performance, robustness against spurious 
instrumental/environmental transients (glitches). 

Instead of using independent estimates of the arrival times at each detector, our DOA 
estimator uses data from (all) detector pairs to estimate the needed propagation delays. 
In addition, it also provides an effective detection statistic, combining the data from 
all detectors in the network, at a remarkably light computational cost. 
DOA reconstruction from arrival-time delay estimation in a network of sensors is a 
well known problem in the technical Literature on Acoustics and Radar (see, e.g., [30] 
for a broad review) . The standard method for time-delay estimation in Gaussian noise 
is (generalized) cross-correlation |31j . which is known to perform reasonably well for 
relatively large signal to noise ratios [33] . Remarkably, the correlation-based estimator 
offers worse performances compared to the XWV in the present context, as shown in 
Sect. 5.2 

This paper is organized as follows. In Sect. 2 we introduce the XWV transform, 
and recall its time-shift properties, which are illustrated for the simplest case of sine- 
Gaussian (henceforth SG) GWBs. In the same section we recall the relationship 
between arrival time delays and DOA. In Sect. 3, we illustrate the proposed XWV 
transform based DOA reconstruction algorithm. In Sect. 4 we discuss the effect 
of noise in the data, and the related DOA reconstruction uncertainties. In Sect. 5 
we present the results of extensive numerical simulations, aimed at characterizing 
the performance of our XWV based algorithm both as a detector and as a DOA 
estimator. The simulations are based on SG-GWBs, but the case of more realistic 
waveforms (including Dimmelmeier and binary merger waveforms) is also discussed. 
In Sect. 6 we include a short discussion of the robustness of the proposed algorithm 
against instrumental/environmental transients (glitches). Conclusions follow under 
Sect. 7. 

2. Rationale. From Cross- Wigner-Ville Transforms to DOAs 

In this section we recall a relevant property of the XWV transform, and illustrate it 
using ideal (sine-Gaussian) waveforms. We further recall the relationship between the 
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arrival time delays and the DOA, for a 3-detectors network, with special reference to 
the LIGO- Virgo Observatory. 

2.1. Cross- Wigner-Ville Transforms 

The XWV transform built from two (analytic, complex) signals ii,2 is given by [34] . 



WMt, f)= J de il \t - - j i2 [t + - j exp {~2Tiif6) (1) 

where * denotes complex conjugation. We recall that the so called analytic signal 
corresponding to a generic real- valued waveform x{t) is 

x{t)^x{t) + iH[x\{t), (2) 

where 'H[x\{t) is the Hilbert transform For xi{t) — X2{t) = x{t)^ eq. ([I]) reduces 
to the well known Wigner-Ville transform of x{t). 



2.2. Time- Shift Property of Cross- Wigner-Ville Transform 
Let To the time-shift operator, such that 

Tg[x]^x{t-e). 

The following property of the XWV transform is easily proved: 



(3) 



■(4) 



Hence, if xi and X2 are the same waveform x{t), except for having different amplitudes, 
and different time-shift (delays), one accordingly has: 



WTgJxi],TB^ [X2]itj) 



= c 



t- 



(5) 



where C is an irrelevant (positive) constant. 



2.2.1. Sine Caussian GWBs To illustrate the practical significance of eq. ([s]) in the 
context of GW detection of unmodeled transients using a network of interferometers, 
we shall refer here to SG waveforms, which have been widespreadly used to model 
GWBs. More realistic transient waveforms will be considered in Section 5.3. Consider 
two SG waveforms, with common carrier frequency /q, time spread T, and initial 
phases </)o, peaked at ii,2, with amplitudes Ai^2, respectively, viz.: 

a;,(t) =^,cos[27r/o(i-tO + '^o]exp[-(t-t,)VT2], z = l,2 (6) 

Under the assumption f^T ^ 1, the analytic counterparts of ([6]) are asymptotically 
given by 

x,{t)^ A,e^^[2TTif^{t-ti)+iMeM-{t-UflT\ i^l,2 (7) 

The XWV spectrum, eq. ([!]) between xi and X2 can be computed in closed form, 
yielding 



W^2{tJ)^W2i{tJ) = (27r)i/2AiA2Texp[-2^2j,2(^_^^)2] 
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• exp 



2^2 



t- 



tl+t2 



It is seen that |Wi2(i, /)| is peaked at 

tl+t2 



t = 



-J / ~ /o- 



exp [-2i7r/(ti - t2)] 



(8) 



(9) 



2.2.2. Realistic Waveforms The above peak localization property of the XWV holds 
true not only for SG waveforms, but essentially for all waveforms modeled by 
oscillatory transients with unimodal envelope, provided the product between the 
(instantaneous) carrier frequency and the envelope duration is a large number. Under 
this respect, the SG waveform is a kind of (worst) limiting case in view of its minimal 
spread property in the time frequency plane. Indeed, the localization property can be 
more marked for other transient waveforms, like, e.g., those numerically generated for 
supernovas or mergers, as discussed in Section 5.3. 



2.3. XWV Spectra, Delays and DO As 

Let us confine for simplicity to the relevant case of the LIGO- Virgo network, which 
consists of the three large-baseline detectors located at Livingston LA (USA) , Hanford 
WA (USA) and Cascina (Italy), henceforth denoted as LI, HI, and V, and labeled by 
the suffix i — 1,2, 3, respectively. In the presence of a GWB, in view of eq. the 
three XWV spectra computed from the data gathered by the LIGO- Virgo network 
interferometers will be (scaled) replicas of the Wigner-Ville transform of the observed 
GWB, exhibiting magnitude peaks at[|] 

t = T,, = ^^^, /-/., =/o, {*,J} = {1,2},{1,3},{2,3}, (10) 

where is the GWB arrival time at dctcctor-i. Knowledge of the three T^- from the 
corresponding XWV peaks allows to retrieve in principle two independent arrival-time 
delays, e.g., 

^13 = n - T3 = 2{Ti2 - T23), t23 = r2 - Ta = 2(^12 - T13), (11) 

from which the DOA, and hence the source location on the celestial sphere can be 
uniquely inferred, as shown in the next subsection. 



2.3.1. DOA from Delays The DOA is easily retrieved from the arrival-time delays 
using the reference system sketched in Figure 1, whose origin is the circumcenter O of 
the triangle whose vertexes are: (1) LIGO-Livingston (LI), (2) LIGO Hanford (HI), 
and (3) and Virgo (V), and whose a;-axis goes, e.g., through LI. In this reference 
system the three detectors have spherical polar coordinates 

{d,^n/2,^ = ipi), i = 1,2,3 (12) 

and are located at: 

ri ~ R{ux cos (fii + Uy sin ipi) , i = 1, 2, 3 (13) 



f We implicitly assume the interferometers' transfer functions as being frequency independent 
throughout the useful band of the sought signals. 
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where = by construction, and R is the radius of the circumference through LI, 
HI and V. Let the source polar coordinates and vector position be ■(? = , (/? = 
and 

f— p(siiii!) cos tpux + sin'd simpuy + cosiJu^), (14) 

respectively, where p is the distance of the source from O. Under the obvious 
assumption where p ^ R, one has 

1^- nl-^ p- Rsin-dscosiips - ^fii) ,i = 1, 2,3. (15) 

whence the delays between the wavefront arrival times at the detectors are 

U] = Ti- Tj = c~^Rsin-ds [cos{ips - ipj) ~ co%{ips ~ , 

{»,J} = {1,2},{1,3},{2,3}, ^'^^ 

where c is the speed of light in vacuum. From the ratio ^ — tiz/t2z one may accordingly 
retrieve ips as follows, 

(1 — ^) cos <(!33 + C COS Lp2 — COS Lpi 



LPs — — tan 



(1 - ^) sin (^3 + ^ sin <^2 - sin ipi 



(17) 



Once ifs has been computed, it can be used in (any of) eqs. (16) to retrieve dg. Note 



that the delays (16 1 do not change upon letting -ds — > it — -ds, yielding the source 
mirror image w.r.t. to the detectors' plane. The above mirror-image ambiguity in a 
3-detectors network is well knowrj|] [5] . 

24. XWV Spectra of Noise 

As a preparation for the next sections, it is important to characterize the key features 
of the XWV spectrum of independent, pure stationary Gaussian noise streams (the 
effect of instrumental transients, aka glitches, will be discussed in Sect. [6|. In this 
case, the time-frequency levels in the XWV spectrum will be random, and their sta- 
tistical distribution, in view of the assumed noise stationarity, will be the same for all 
(discrete) times. 

The first two moments of the above distribution can be computed analytically with 
relative ease ,37 . In particular, for all (discrete) frequencies the average value is zero, 
and the variance exhibits a piecewise linear dependence on frequency, as sketched 
in Figure 2. The maximum variance occurs at / = /s/2, where fg is the sampling 
frequency, and its value depends on the details of the XWV implementation (size, 
windowing), and the noise level in the data streams (see Appendix for details). It is 
thus expedient to equalize the XWV time-frequency levels, so as to obtain a uniform 
(flat) XWV spectrum for pure-noise data streams. To this end, we merely scale the 
XWV level in each time frequency pixel to the (computed) standard deviation of the 
XWV level in that pixel. 



3. Estimating DOAs from Discrete XWV Spectra of Noisy Data 

In practice, the XWV spectra will be computed in discrete form |34j. yielding two- 
dimensional (complex) arrays, rather than continuous time-frequency functions over 

§ The mirror-image ambiguity can be resolved, in principle, from knowledge of the detectors' pattern 
functions, featuring different responses in the i9 = i9s and i9 = tt — i?s directions. 
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To minimize the effect of time-discretization error it is convenient to estimate the 
independent delays corresponding to the largest available baselines, i.e., in our case, 
ti3 (Ll-V) and t23 (Hl-V). 

Also, in the presence of noise eqs. (11) used in (16) will provide a mere estimate of 
the DOA, whose goodness will basically depend on the available signal to noise ratio, 
which affects the accuracy whereby the XWV peaks can be identified. 
A simple algorithm for seeking peaks in the three LIGO- Virgo XWV spectra which 
are consistent with the constraints 

= 2|ri2-Ti3| <c-i|r23| 

= 2|r23~Ti2| <c-i|ri3| (18) 
= 2|r23-Ti3| <c-i|ri2| 

expressing the obvious requirements that the wavefront propagation delay between 
two detectors cannot exceed the limiting value corresponding to propagation along 
the line-of-sight direction between the detectors, can be now formulated. The algo- 
rithm uses the three (discrete, noisy) XWV spectra to construct a grid in the time 
delay plane (^13,^23), and assign different levels R to its nodes: 




initialize all time-delay grid node levels to zero 
for all time-frequency pairs (ri2,/i2) in W12 

for all time-frequency pairs (Ti3,/i3) in W13 such that: 

2|ri2 - T13I < c"^|f23| and /13 = /12 

for all time-frequency pairs (723, 723) in W23 such that: 
2|Ti2 ~ 723I < """1^131 and /23 = /12 

accumulate level i? = i?+|VKi2(ri2, /i2)VKi3(ri3, /i3)W^23(r23, /23)| 
at grid node {ti3 = 2(ri2 - 123), i23 = 2(Ti2 - T13)} 

end for 

end for 
end for. 



A candidate direction of arrival is obtained by taking the highest-level grid-node in 
the (^13,^23) plane subset defined by the further constrainijji] 

1^12! = 1^13-^231 < C-M^^12| (19) 

The highest level in the grid can be used both as an estimator of the DOA, and as a 
detection statistic, whose performances will be discussed in Sect. 5.1. 
Note that the proposed algorithm is coherent, in the sense that it produces a single 
detection statistic by combining the data from all detectors in the network. It also 
inherits the typical features of coincident tests: the outermost loop enforces frequency- 
consistency, while the two inner loops enforce time-delay admissibility. 
Note also that the XWV spectra will display sensible peaks only if the waveform 
gathered by the different detectors are consistent in shape. This suggests that the 
algorithm will be robust against (independent) instrumental disturbances, as further 
illustrated in Sect. 6. 



As shown in Sect. 4, the bound in il9l can be made sHghtly tighter. 
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3.1. LIGO-Virgo Network Directional Response under XWV Based Algorithm 

In the absence of noise, the above algorithm will produce a peak in the time-delay grid 
whenever a GWB is observed by the LIGO Virgo network. The peak will be located at 
a node whose time-delay coordinates correspond to the DOA {ds, <Ps)- This peak will 
be well localized provided the duration of the transient signal is substantially shorter 
than the minimum graviton flight-time between detectors. 

In view of the bilinear nature of the XWV, it can be argued that the peak height will 
be proportional to the squared product of the three detectors' pattern functions along 
that direction. This quantity, normalized to its maximum, and denoted henceforth 
as ^{^, f) describes the directional response of the proposed GWB detector/DOA 
estimator, and is plotted in Figure 3 for the LIGO-Virgo network, for circularly 
polarized GWs. We checked numerically that the expected (normalized) levels of 
the time-delay grid peaks, reproduce those computed from the function $ for each 
DOA in a iy9 grid of 50 x 100 points (using 10"^ noise realization for each DOA). 
The quantity: 



47r 47r 



smMMip (20) 

*(t5,I/3)>'I'mi„ 



expresses the fraction of the (unit) celestial sphere where the (normalized) directional 
response of the proposed detector/estimator exceeds the threshold value ^mvm and is 
displayed in Figure 4. 

It is seen, e.g., that roughly 50% of the celestial sphere is covered with <I>,„i„ > .2 by 
the LIGO-Virgo network, using the proposed algorithm. 

4. DOA Reconstruction Uncertainties 

Uncertainties in the DOA reconstruction stem from a twofold origin: the discreteness 
of the time-delay grid, due to the discrete implementation of the XWV spectra (finite 
time resolution), and the additive noise in the data (see discussion in Sect. 5.2). 
In order to translate the effect of systematic and statistical errors in the estimated 
delays into uncertainty ranges in the estimated DOAs, it is expedient to introduce the 
projection which maps the DOA polar angles (i?s, ips) into a point {xg, Vs) of the disc 
(with center O and radius R) going through the detector, viz[^ 

Xs = i? sin COS (^s , — i? sin i?^ sin (p^ . (21) 

The formula which relates the {xs, Us) projection to the arrival-time delays is obtained 



from eq. ( 16 ) 



Vs 



(ti3 - ^23)sinv73 - ti3sin(^2 
sin((/72 - fs) + smip2 - sin 993 

(22) 

^13 COS (p2 - {ti3 - ^23) cos (p3 - ^23 

- sin((y52 - "^3) + sin (^2 - sin (/33 

Equation (22 1 is a Zmear transformation, relating not only the coordinates {xs,ys) to 
the delays (ii3,i23), but also the uncertainties Sxg, Si/s to the delay errors Sti^ and 
6t23. Thus, under the simplest assumption where these latter are independent and 
identically distributed, the uncertainty region in the (ii3,t23) plane is a circle, and 

^ We recall that the x-axis goes through detector-1 (LIGO-Livingston). 
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the corresponding uncertainty region in the (xs,ys) plane is an ellipse. Notably, the 
shape of this latter is translation-invariant across the circle + < R^, i.e., DOA 
independent. 

The ratio between the uncertainty areas in the {xs,ys) and (ti3,i23) planes is given 
by the Jacobian of the transformation (22 1, viz.: 

J = -. ""^^-^ ' , (23) 

sin(/33 - smLp2 + sm((p2 - fs) 

which is also DOA-independent. 

By back-projecting the uncertainty ellipse onto the sphere of radius R centered at O, 
we obtain a DOA-dependent uncertainty region. This is illustrated in Figure 5, for a 
few representative cases. The ratio between the area of the uncertainty region on the 
celestial sphere, and the area of the uncertainty ellipse in the {xs, j/s) plane is displayed 
in Figures 6a and 6b as a function of cps, for various values of d^. For i?s ~ 0, this 
ratio is close to unity, whatever ips- On the other hand as 1?^ '^/'^i the ratio blows 
up, and its dependendency on ip^ becomes more and more evident. Such a behaviour 
had been already noted in, e.g., [TT], |13) . 



5. Numerical Experiments 

In order to check the performance of the proposed algorithm, we run a series of 
Monte Carlo simulations. The simulations use time-discretized GWBs and glitches 
injected into white (independent) random Gaussian sequences, to represent the three 
interferometer data. In the case of GWBs, the delays are chosen according to the 
assumed source location. 

Our XWV engine uses data chunks 2048 time-samples wide to produce a 1024 x 1024 
time-frequency nodes XWV transform, using Pei-Yang fast algorithm [38] . The 
sampling frequency is 4 KHz. The data are accordingly decomposed into half- 
overlapping chunks 2048 time-samples wide (we use a plain rectangular windowing 
function), in order to use a fixed number of time samples to compute each time- 
frequency samples. The resulting discrete XWV spectrum spans the time range 
between samples #513 and #1536, and the frequency range between and 1000 7720 
As already mentioned, the XWV values are equalized so that in the absence of signals 
their first and second moment are and 1, respectively. 

Now, even in the absence of a signal, the levels produced by our algorithm in the time- 
delay plane (ti3,t23) grid-nodes will be non-uniform, due to the different number of 
(noisy) time-frequency XWV values mapped into each node. The average and standard 
deviation in the (ti3,t23) plane for pure-noise (stationary, Gaussian) data are shown 
in Figure 7. We accordingly equalize the levels, by subtracting the above average, and 
dividing the result by the above standard deviation, so that in the absence of signals, 
the grid-node levels in the time-delay plane will have zero average and unit variance. 
For each injected waveform, we generated 10** different noise realizations, to test the 
statistical properties of the proposed algorithm, both as a detector and as a DOA 
estimator. The waveforms were parameterized by their intrinsic signal to noise ratio 
(SNR), defined by 

. ^Kss ^ {I[hl{t) + hl{t)]dt}'^' 

^ N N ' ^ ' 

+ Note that when using analytic signals for co mpu ting discrete versions of the XWVT, the minimum 
sampling rate must be twice the Shannon rate |34| 
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being the (two-sided) power spectral density of the stationary white(ned) Gaussian 
noise component, assumed for simplicity the same in all detectors (the effect of glitches 
will be discussed in Sect. 6). 

The results of our simulations are summarized below. 

5.1. Detection Performance 

The performance of our algorithm as a detector are illustrated in Figures 8 and 9. The 
detection statistic is the level of the highest peak in the time-delay grid. Figure 8a 
and 8b display the false alarm (continuous line) and false dismissal probabilities (the 
dashed lines, corresponding to different values of the intrinsic SNR (Sh) as functions 
of the detection threshold 7, for DOAs corresponding to the maximum = 0.705rad, 
If = 5.073rad) and the minimum {•& — 0.800rad, tp — l.lOOrad) of the network pattern 
function in Fig. 3. Figures 9a and 9b show the receiver operating characteristics, i.e., 
the detection probability vs the false alarm probability, for fixed values of the intrinsic 
signal to noise ratio, Sh, for a DOA corresponding to the maximum of the network 
pattern function in Fig. 3. 

5.2. DOA Estimation Performance 

As already mentioned, the finite time resolution implies that the estimated delays 
are affected by a systematic uncertainty which can be twice the XWV time-step 
St. The noise in the data entails that estimated delays spread around the actual 
delays in a signal-to-noise dependent way. This is illustrated in Figures 10 and 11. 
The estimate is always unbiased, whenever the signals are shorter than the minimum 
graviton flight time between the detectors. Figure 10 displays the standard deviation 
of the estimated delays (average between the two) as a function of the intrinsic SNR, 
for DOAs corresponding to the maximum = 0.705rad, = 5.073rad) and the 
minimum (t? = 0.800rad, (p = l.lOOrad) of the network pattern function in Fig. 3. 
In a log-log scale, both curves show the same slope, corresponding to an exponent 
« —1.4. Figure 11a displays the empirical distribution of the estimated delays for 10"* 
different noise realizations, for a DOA corresponding to the maximum of the network 
pattern function in Fig. 3, for two different values of the intrinsic signal to noise ratio. 
It is interesting to compare Fig. 11a to Fig. lib, where a standard correlation-based 
time-delay estimator [33] has been used to retrieve the two propagation delays. Our 
XWV-based estimator is seen to offer distinctly better performances. 

5.3. Realistic Waveforms 

As anticipated in Section 2, the XWV transform peak localization property holds not 
only for SG waveforms, but for general transient waveforms. This is further illustrated 
in Figures 12 to 15. 

Figure 12 (top) shows two copies of a typical supernova GWBs, belonging to the family 
computed by Dimmelmaier and co-workers |40| . with a time shift of 82 time samples 
(corresponding to 20.5 ms, at our sampling frequency), together with their XWV 
transform. The XWV is identical to the Wigner transform of the GWB waveform, 
except for the time-shift given by eq. ([5|. Accordingly, the XWV peak in Fig. 12 (bot- 
tom) is localized at the midpoint between the peak times of the two waveforms in Fig. 
12 (top). Figure 13 (left) shows the time-delay grid histogram when this waveform is 
emitted by a source located in the direction of maximum network sensititivity in Fig. 
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3, for Sh — 100. It can be seen that, not unexpectedly, the locahzation properties in 
the delay plane are even better than for Gaussian waveforms, in view of the larger 
time-bandwith figure of the Dimmelmaier waveform. 

Figures 14 and 15 are similar to Figs. 12 and 13, except that the waveform here is 
that of a typical binary merger [41] . 

6. Glitch Rejection 

By construction the proposed detection/localization algorithm should be robust 
against spurious transients of environmental/instrumental origin (glitches). 
We may expect that glitch induced false detection may occur only in the (unlikely) 
case where each detector shows a glitch in the analysis window, such that the mutual 
delays are consistent with an acceptable DOA, and the (indipendent) glitch waveforms 
are consistent in shape. 

In order to illustrate these features, we consider first the no-GWB case where a glitch 
occurs in the data of each of the three interferometers, the three glitches being differ- 
ent, but with delays consistent with an admissible DOA. 

To this end, we used a set of = 7 visually different waveforms, shown in Figure 
16 (top), from the catalogue of "typical" LIGO glitches compiled by P. Saulson [42] . 
All glitches in the set were scaled to unit norm, and time-shifted so as to bring their 
envelope peaks to coincidence. The (normalized) pairwise correlation coefficient of the 
selected glitches, which provides some quantitative measure of their (dis)-similarity, 
does not exceed 0.62, with mean and median values of 0.195 and 0.105, respectively. 
The correlation coefficient histogram is shown in Figure 16 (bottom left). 
From the above glitch set, we formed (all) 35 triplets of different waveforms, and 
computed the related X-Wigner transforms and the peak-levels in the time-delay grid 
produced by our algorithm. For each of glitch-triplet (51,52,53) we also computed the 
geometric mean of the time-delay grid peak-levels for the three cases where the data 
from all interferometers contain the same waveform 5^ , i = 1,2,3. This quantity 
was used to re-scale the time-delay grid peak- level for that glitch-triplet. 
The histogram of the rescaled time-delay grid peak-levels for the 35 different glitch 
triplets considered is shown in Figure 16 (bottom right). The largest (scaled) peak 
level was 0.51, with a median value of 0.12 and a mean of 0.17. 

These, admittedly limited, results illustrate the waveform consistency test capabilities 
of the proposed algorithm. 

We next consider the case where GWBs and glitches co-exist in the data. In these 
further simulations we used SG glitches and GWBs, with glitch parameters (center fre- 
quency, peak position and carrier frequency) generated randomly and independently 
in each detector. The pertinent results are illustrated in Figures 17 to 20. These 
figures show the noisy waveforms (left column), the density maps of the XWV trans- 
forms (mid column), and the (normalized) level map in the (^12,^13) time-delay grid. 
In Fig. 17 we consider the simplest case where the GWB data in a single detector (HI 
in this case) are corrupted by a single glitch in the analysis window. The GWB signal 
to noise ratios are 22.16 (LI), 23.13 (HI) and 31.85 (V), for a (circularly polarized) 

* These correspond to a source equidistant from all detectors radiating the waveform gi, with all 
detectors exhibiting the same response in the source direction. The last assumption is unrealistic, 
but is irrelevant for the present purpose). 
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source with (5/^ = 50 at i? = 2.58rad and (f = 2.71rad. The glitch signal to noise ratio 
in HI is 22.59. The ghtch shows up clearly in the HI data, and produces evident arti- 
facts in the Hl-V and Hl-Ll XWV spectra. Nonetheless, its effect on the time-delay 
level map is almost negligible. 

In Fig. 18 wc consider the (unlikely) case where the data in each detector arc cor- 
rupted by (single) glitches in the analysis window, with SNR values of 11.5 (LI), 12.00 
(HI) and 16.53 (V). None of the glitches has a significant overlap with the GWBs; 
nonetheless, they produce artifacts in all XWV transforms. Also in this case the effect 
of these artifacts on the detection/localization properties is negligible, as seen from 
the time-delay level map. 

Not unexpectedly, the localization performance deteriorates significantly in the rather 
extreme situation where glitches overlap the GWBs in the data. When this happens 
in such a way that true peaks in the XWV transforms are no longer resolvable from 
spurious ones, the time delay level map topography is substantially blurred in the 
neighbourhood of the true delays, resulting into more or less severe localization errors. 
This is illustrated in Figs. 19 and 20. 

In Fig. 19 wc have a (single) GWB-ovcrlapping glitch in HI with SNR=16.28. In Fig. 
20 each detector is affected by a GWB-overlapping glitch. The glitch SNR values in 
Figs. 19 and 20 are the same as those in Figs. 17 and 18. A sensible distortion in the 
XWV transforms is observed, entailing a sensible error in the estimated delays. 

7. Conclusions 

We presented a simple, computationally light and fast algorithm for detecting 
short unmodeled GWB in a network of three interferometric GW detectors, and 
estimating the related DOA, based on XWV spectra. The algorithm is reasonably 
performant, and nicely robust against spurious transients (glitches) of instrumental 
origin corrupting the (otherwise Gaussian) detectors noise floor. 

It does not provide waveform reconstruction; this latter, however can be accomplished 
in principle off-line, once the DOA has been estimated. 

Generalization to larger networks, and other potentially interesting waveforms (e.g., 
chirps) is relatively straightforward. Such extensions will be explored in a forthcoming 
paper. 

Based on the above preliminary results, we suggest that the proposed algorithm may 

be used as a quick-and-(not-so)-dirty on-line data sieving tool. 

A quantitative comparison with existing GWB detection/DOA estimation algorithms 
in terms of efficiency and computational burden will be the subject of future 
investigation. 
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Appendix - XWV Moments. 

Let 

L 

Wf^g{n,m) = f{n + k)g*{n- k)exp[-'iiN-\mk]. (25) 

k=-L 

the discrete version of tlie XWV spectrum, where n and k are the discrete time and 
frequency index. Note that, formally, 

WfJno,m.)^DFT,i2m) (26) 

where x = f{no + k)g*{no — fc). This shows that while the time index n spans the 

range {—L, L) n N, the frequency index m spans the range (— L/2, L/2) n N. 

Let 

i> = iy + iHiy (27) 

the analytic version of the background noise, and denote as i'(fc) and ^^{k) the (real- 
valued) samples, of the noise and its Hilbert transfom. 

For zero average Gaussian white noise with (two sided) power spectral density Wq^ 
the spectral power density of analytic noise is 

= WoC/[^(m)], (28) 

where U{-) is Heaviside's function and 9 — 2Tim/N £ (— 7r,7r), = 2_L + 1 being the 
number of DFT frequency samples. 

It is a simple task to show that the first moment of the XWV is zero, in view of the 
assumed independence of the (Gaussian) noises in different detectors. 
We now compute the second moment, viz.: 

L L 

= E E ^K(" + fcK("-fc)^:(n + p)/ia(n-p)]e-^2(fe-rte(™)^ 

k=-Lp=-L 
L L 

= E E ^K("+fc)^:("+p)]^K(n-p)M:(n-fc)]e-^2(fc-p)eM 

k—~L P——L 

= E E RL^Ak-p)e'''^'~^^'^"'^ (29) 

k——L P——L 

where Va and fj,a are built from independent, zero average, white(ned) Gaussian noises 
pertinent to different detectors, but assumed as having the same power spectral density 
Wo, and Ri,^^^^{h — k) = E[va{h)i'l{k)\ denotes the autocorrelation, function of 
analytic noise. The inner summation in ( |29[ ) can be extended to cx), and the Wiener 
Khinchin theorem can be invoked to prove that, 

L oo 

a2(n,™)« Y E <,.„(fc-p)e-'^('=-^)^(") = 

k=~L p=~oo 
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k=-L 

where m e [-L/2, L/2] n N. 



TT 



32W^\m\ 
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(30) 
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Figure 2 - Normalized variance of XWV of stationary Gaussian noise 
vs. frequency (fixed time-slice). 
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Figure 5 - Uncertainty areas on the celestial sphere 
and the detectors' plane for some source positions. 
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Figure 6a - Ratio between the uncertainty areas on the celestial sphere (AST) 
and the detectors' plane (AQq) vs. source position (&,(p). 
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Figure 7 - Expected value (left) and standard deviation (right) of 
proposed network detection statistic. Gaussian noise only. 
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Figure 8a - False alaim and false dismissal probability vs. threshold, for various values 
of the intrinsic SNR (5^ in eq. 24). Source at maxima of the pattern function in Fig. 3. 
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Figure 8b - False alarm and false dismissal probability vs. threshold, for various values 
of the intrinsic SNR (Sjj in eq. 24). Source at minima of the pattern function in Fig. 3. 
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Figure 9b~ROCs, for various values of the intrinsic SNR (B,, in eq. 24). 
Source at maxima of the pattern function in Fig. 3. 
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Figure 10 - Standard deviation (in time bins) of estimated arrival times 
vs. intrinsic SNR (5y, in eq. 24). Source at maxima (dashed line) 
and minima (solid line) of the pattern function in Fig. 3. 
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Figure 1 1 a - Nomialized binning (histogram) of estimated time delays, 
for two representative values of the intrinsic SNR: 5/, =15 (left) and Sy, =35 (right). 
Source at maximum of pattern function in Fig. 3. 




Figure 1 lb - Same as Fig. 1 la, but a standard correlation- based 
time delay estimator has been used to retrieve the propagation delays ?y,, / 
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Figure 13 - Normalized density (left) and surface plot (right) of time-delay 
grid levels. Supernova core-collapse GWB (Dimmelmaier). 
Source at maximum of pattern function in Fig. 3. 
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Figure 15 - Normalized density (left) and surface plot (right) of time-delay 
grid levels. Binary merger GWB (Baker). 
Source at maximum of pattern function in Fig. 3. 
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Figure 1 6 - Top: set of typical glitches (amplitude vs time bin) from [42], scaled to unit 
energy. Bottom left: histogram of related pairwise correlation coefficients. Bottom right: 
histogram of normalized time-delay grid-peak levels for triplets of different glitches. 




Figure 17 - GWB and single glitch (non overlapping) in one detector (HI). Left: 
waveforms; Mid: XWV spectra; Right: normalized density plot of time delay grid levels. 
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Figure 18 - GWB and single glitch (non overlapping) in each detector Left: wavefijrms; 
Mid: XWV spectra; Right: nonnalized density plot of time delay grid levels. 
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Figure 20 - GWB and single glitch (overlapping) in each detector. Left: wavefijrnis; 
Mid: XWV spectra; Right: normalized density plot of time delay grid levels. 



